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A path- integral representation is constructed for the Jahn- Teller polaron (JTP). It leads to a pertur- 
bation series that can be summed exactly by the diagrammatic Quantum Monte Carlo technique. 
The ground-state energy, effective mass, spectrum and density of states of the three-dimensional 
JTP are calculated with no systematic errors. The band structure of JTP interacting with disper- 
sionless phonons, is found to be similar to that of the Holstein polaron. The mass of JTP increases 
exponentially with the coupling constant. At small phonon frequencies, the spectrum of JTP is flat 
at large momenta, which leads to a strongly distorted density of states with a massive peak at the 
top of the band. 

PACS numbers; 71.38.-f-i, 02.70.Lq 



The renewed interest in the Jahn- Teller polaron (JTP) 
problem was sparked by its possible relevance to the 
colossal magnetoresistance phenomenon in manganese 
oxides [|l|-g. Although cooperative properties of JTPs 
are of most importance, the solution of the single-JTP 
problem is a necessary step towards the comprehensive 
understanding of the physics of these and other materi- 
als with Jahn- Teller ions. JTP also bears a significant 
fundamental interest and it has been a subject of in- 
tense theoretical research throughout decades, see, e.g., 
P-p[. However, the full quantum-mechanical treatment 
of this problem is difficult, and no exact solution has 
been found. Even the single-site JTP does not allow an- 
alytical solution. The latter fact has so far prevented 
the development of a strong-coupling perturbation the- 
ory, so successful in the theory of the Holstein polaron 
M. Numerical methods, such as exact diagonalization 
and density-matrix renormalization group, face the po- 
tential difficulty in dealing with a very large Hilbert space 
of the two phonon modes which are the essential feature 
of the JTP. Quantum Monte Carlo studies reported so 
far have treated the phonons only classically ||l^ . 

In this paper, a method is presented that allows the 
full quantum mechanical solution of the JTP problem. 
The method combines algorithms developed in the Monte 
Carlo studies of the Holstein ||ll|,|2| and Frohlich [|3) po- 
larons. It is based on a path-integral representation of 
JTP partition function, analytical integration over the 
phonon coordinates |14| , |ll| , and numerical integration 
over the electron coordinates. The latter is done exactly, 
with the Diagrammatic Quantum Monte Carlo technique 
[ jl3[ . The use of free boundary conditions in imaginary 
time |1^ allows calculation of the polaron effective mass 
and the entire spectrum. The method is formulated in 
continuous time and no systematic errors are introduced. 
As a result, the ground-state energy, mass, spectrum, and 
density of states of JTP are calculated exactly on an in- 



finite lattice in the most difficult three-dimensional case. 
We found that the properties of JTP are very similar to 
that of the Holstein polaron (HP). JTP is as heavy as HP. 
It also displays some anomalous properties, again similar 
to HP. In particular, in the case of dispersionless phonons 
and low phonon frequency, JTP has a flat spectrum at 
large momenta and a very distorted density of states. 

The simplest JTP is an electron that hops between 
doubly degenerate eg levels and interacts locally with a 
doubly degenerate phonon mode. The model Hamilto- 
nian reads H 
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Here cjj^j creates an electron on site n on orbital a {a = 1 
or 2 is the orbital index), t is the intersite hopping ma- 
trix element, (nn') denotes pairs of nearest neighbors; x^ 
and y„ are the two phonon displacements on the nth site, 
Pxn,yn — ~*^^a;n-yni ^^ ^^ ^hc ioulc mass; and n is the in- 
teraction parameter with dimensionality of force. The 
energy of the atomic level is set to zero. The model is 
parameterized by the phonon frequency uj and by the di- 
mensionless coupling constant A = k^ /{2Muj'^zt) where 
z is the number of nearest neighbors. For the simple 
cubic lattice, considered in this work, 2 = 6. Two sim- 
plifications are introduced in the model (E])-(0). First, 
the electron hopping is isotropic, and is diagonal in the 
orbital index. The latter leads to any change of orbital 
index taking place only when an x-phonon is emitted or 



absorbed. Averaging over the phonons eliminates all the 
amplitudes with any odd total number of emission and 
absorption acts. As a result, the off-diagonal in orbital 
index elements of the polaron density matrix are iden- 
tically zero. One can show that such matrix elements 
determine the splitting between the two polaron bands. 
The conclusion is that in the model (l])-(0), the two po- 
laron bands are degenerate in the whole Brillouin zone, 
and there exists only one ground-state dispersion rela- 
tion Ep . The second simplification is the absence of the 
phonon dispersion. It is assumed that each site is sur- 
rounded by its own vibrating cluster and the neighboring 
clusters share no common ions. Such a choice was made 
to bring the model close to the Holstein model iQ. The 
difference between the two models lies exclusively in the 
interaction term (Q). This fact will enable us to make a 
reasonable comparison of JTP and HP. One should add 
that the full Hamiltonian (||) is invariant under the con- 
tinuous transformation: rotation in the plane (x, y) by an 
angle 4> (the same for all sites) and simultaneous rotation 
in the plane (ci,C2) by (/'/2. The invariance implies the 
equivalence of the x and y modes in the sense that the 
elastic and interaction energies associated with the two 
modes must be equal. The number of excited phonons 
are equal too. These general results serve as independent 
checks of numerical data. 

Central in the proposed method is the path-integral 
representation for the density matrix p of the model (0)- 
(0). For a small imaginary-time interval Ar — /3/L, 
where j3 = [ksT)^^ is the inverse temperature and 
L ^ 1, one obtains, up to the first order in Ar: 

p(Ar) = (r', a'; {x'^}, {y^jle^^^^lr, a; {x„}, {y^}) 

= [5rr' 5 aa' + ^T ndrr' 5aa' Xr 

-HArMaa' Y^ <5r,r'+„-„']e^?'-\ (5) 
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^ph being the phonon action. Here a — 1, (2) when 
a — 2,(1), i.e. a is 'not' a. The full density matrix 
is obtained by multiplying p(At) by itself L times, in- 
troducing integration over the internal coordinates, and 
taking the L — > oo limit. In doing so, we do not impose 
the periodic boundary conditions in imaginary time but 
rather leave initial r(0), a(0) and final r(/3), a(/3) coordi- 
nates of the electron arbitrary. [Except that a(/3) = a(0) 
always holds.] The boundary conditions on the phonon 
coordinates are twisted, XniP) — a;n-r(/3)+r(o)(0), and 
the same for y, as described in ||l5]. The next step is to 



integrate over Xn{T) and yn{T). This is a Gaussian inte- 
gration and it can be done analytically. Here one faces 
an important difference between the two phonon modes. 
Displacements ynir) interact with the electron density 
and do not change the electron coordinates. Therefore, 
for the oscillators j/n, the electron is simply a source of 
external force, the value and sign of which depend on the 
position and orbital index of the electron, cf. the sec- 
ond term in Eq. (^). The problem reduces to uncoupled 
oscillators in a time-dependent external field. The inte- 
gration over ynir) can be done by Feynman's methods 
p3 to yield the factor e"^" in the density matrix, where 
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is the retarded action induced by the y mode. The second 
term in the function G{t — t') involves the net shift of the 
electron path Ar = r(/3) — r(0) and the sign of the time 
difference sign(r — r') = ±1, |13|. The simple form (0) 
of G is valid in the limit e^^'^ ^ 1 that is assumed here- 
after. It is important that Ay is an explicit functional 
of the electron path and can be calculated straightfor- 
wardly once the functions r(T) and a{T) are specified. 
Note, that the y-mode 'favors' the same orbital index 
throughout the whole electron path and 'dislikes' orbital 
changes. 

Unlike the y-mode, the x-mode changes the orbital 
index of the electron. It therefore acts like the kinetic 
energy, but in the 'orbital space', compare the last two 
terms in brackets in Eq. (|^) . The only complication is that 
the rate of orbital change is not constant but depends on 
the local displacement Xrir). Nonetheless, the x-mode 
and kinetic energy can be treated in the same manner, 
and this is done as follows. Upon the self- multiplication 
of /9(At), a variety of terms with different powers of (At) 
appear. A term with n site changes ('kinks') and m or- 
bital changes has the weight 
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Here r^ is the time of the sth orbital change (s — 
1, . . . , to), v{ts) is the electron position at this time, and 
Xj.t^\{Ts) is the x-displacement at this site at this time. 
Integration over Xn(T) in (pi) can now be performed by 
introducing fictitious sources, calculating the generating 
functional and differentiating it m times. For any odd 
TO. the result is zero, since the x action is even under the 



global sign change Xn(T) -^ —Xn{T). (This is the above- 
mentioned averaging that leads to the complete degen- 
eracy of the two polaron bands.) For even m, the time 
moments r^ should be combined in pairs, and each pair 
contributes the factor k^G{ts ~ Ts') with G from Eq. (H). 
Because G is always positive, the a:-mode 'favors' a large 
number of G- lines on the electron path, i.e., a large num- 
ber of orbital changes. Such a tendency is opposite to 
that of the y-mode. The competition between the two 
modes leads to a dynamical balance and to some average 
number of orbital changes per unit imaginary time. Note 
also that the same function (pi) determines both x and 
y contributions to the density matrix, which reflects the 
equivalence of the two modes. 

Now recall that each term (g) must be still integrated 
over all the electron positions and orbital indices, i.e., 
over the time instances of its kinks and orbital changes. 
Such an integration is denoted below as (dr)" and (dr)™. 
In the L — > cxD, At — *■ limit this leads to the path- 
integral representation of p for JTP: 

oo oa .p .fj 
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One cannot proceed analytically any further. However, 
all the integrands are positive-definite, which suggests 
simulation of p by Monte Carlo (MC) methods. The 
obvious difficulty is that an infinite series of integrals 
of ever-increasing dimensionality, rather than just one 
integral, have to be evaluated. Nonetheless the recently 
developed Diagrammatic MC method ||I^ is capable of 
doing such an integration. The method works as follows. 
Each (n, to) term in the sum is represented by a diagram 
(path) with n kinks and m/2 G-lines, see Fig. |l|. The MC 
process browses the configuration space by means of the 
two elementary subprocesses: (i) inserting and removing 
kinks, which changes the dimensionality of integration 
by 1 (n -^ n ± 1); (ii) attaching and removing G-lines, 
which changes the dimensionality by 2 (to, -^ 77i±2). The 
central idea of the method is the balance equation for the 
two subprocesses |13 . Let Nk be the number of kinks of 
a given sort, N^ < n. Next, let R{t') be the normalized 
probability density with which one chooses the position 
for the new kink. For example, one may decide that 
all the times in the [0,/3] interval are equivalent, hence 
R{t') = 1//3 = const. In the reciprocal removing process, 
one may decide that any of the Nk + 1 existing kinks are 
removed with equal probability {Nk + 1)~^ . The resulting 
balance equation reads: 
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FIG. 1. A typical real space - imaginary time diagram 
(path) with n = 8 kinks (5 right and 3 left ones, which makes 
the net shift of the path be Ar = 5 — 3 = 2 lattice sites), 
and m = 12 orbital changes (6 G-lines). The two different 
orbital states are shown as black and white segments of the 
path. Each kink contributes to path's weight the factor t, and 
each pair of orbital changes connected by a dashed line the 
factor k^G{t' — t"). Additionally, there is an overall factor 
e " induced by the y phonon mode. Insertion or removal of 
a kink at time r' shifts the whole segment t' < r < /3 by one 
lattice site. Ar also changes. [Recall that Ar is an explicit 
parameter of G, Eq. (H).] Attaching or removing of a G-line 
at times t',t" changes the orbital index ('color')_of the path 
at r < T < T ' 



This affects the action Ay, Eq. (M). 



from where the acceptance probabilities P follow in the 
usual manner pq ]. The main feature to note is that 
the elements of the phase space associated with the two 
sides of the equation have the same measure: Wnm brings 
(dT)"+™, and R{t') — 1/ (i adds one more (dr) because 
i? is a probability density. The measure of the right- 
hand-side is also (dr)""'"™"'"^. That makes both transition 
probabilities be of the same order, and renders the whole 
process meaningful. In the case of G-lines, let S{t',t") 
be the two-dimensional probability density to attach a 
new G-line at times r' and r". While removing, each of 
Ng + 1 G-lines may again be chosen with equal proba- 
bility. The balance equation reads 
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Again, the measures of the phase space elements are the 
same because now S adds (dr)^ . A possible choice for S 
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FIG. 2. Inverse effective mass of the three-dimensional 
JTP (circles) and HP (squares, adopted from [12]) for 
hu) — 1.0 i, in units of I/ttiq = 210^/%^, ao being the lattice 
constant. Statistical errors are smaller than the symbols. 
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FIG. 3. Density of states of the three-dimensional JTP 
for hu> — 1.0 i and A = 1.3. The polaron bandwidth is 
W = 0.130(5) t. The dotted line shows the free-particle DOS 
with the same W. 



is S{t',t") = {huj/l3)expi~hLu\T' -t"\)). The MC pro- 
cess that follows the rules (|l2|)-(|l3|) generates a Markov 
chain of paths distributed in accordance with Eg. (11). 
On such an ensemble various polaron properties can be 
measured with the standard Metropolis rules fl^ ]. 

One physical quantity that can be calculated with this 
method is the ground state energy of the polaron Eq = 
{-W-^idWnm/dl3)) |ll]jl|. Here we present MC data 
for two other properties: the effective mass and density 
of states (DOS). The mass is calculated as the diffusion 
coefficient of the polaron path, m~^ = {ph'^)~^ {{Ar)1) , 
[ fi2[ . The inverse mass of the three-dimensional JTP is 
shown in Fig. 0. After the initial weak-coupling growth 
m/mo = 1 -|- const • A, a transition to the small polaron 
state takes place at A ~ 1.1 — 1.3, and after that the 
mass increases exponentially with coupling. The com- 
parison with the Holstein polaron (HP) shows that both 
masses behave similarly. The polaron spectrum is cal- 
culated with the formula Ep — Eq = — /3^^ In(cosPAr), 
[ p2[ , and DOS is obtained by integrating Ep over the 
Brillouin zone. The resulting DOS appears to be strongly 
distorted in the adiabatic regime, hco — 1.0 i, see Fig. ||. 
The polaron spectrum is flat in the outer part of the 
Brillouin zone due to hybridization with dispersionless 
phonon modes [|l7| , which results in a massive peak in 
DOS at the top of the band. The van Hove singulari- 
ties are invisible because they are absorbed by the peak. 
All these features are very similar to that of the Holstein 
polaron (l2|. The similarity suggests that such peculiar 
features of the band structure are common to polaron 
models with local interaction and dispersionless phonons 
in the adiabatic regime. 

In conclusion, we have developed a path-integral rep- 



resentation for the Jahn- Teller polaron and used the Di- 
agrammatic Monte Carlo method to simulate the series 
expansion for its density matrix. The ground state en- 
ergy, effective mass, spectrum, and density of states have 
been calculated with no systematic errors. JTP has been 
found to be very similar to the Holstein polaron and, 
therefore, to possess some anomalous properties, in par- 
ticular a strong distortion of the density of states in the 
adiabatic regime, huj < t. Such a distortion may be a 
consequence of the short-range electron phonon interac- 
tion and dispersionless phonon modes. 
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